Hexatic phase in the two-dimensional Gaussian-core model 



O 
(N 

in 



o 



I 

o 
o 



00 
(N 

00 

q 
o 



X 



Santi Prestipino^ *], Franz Saija^ ft], and Paolo V. Giaquinta^ [j] 

^ Universita degh Studi di Messina, Dipartimento di Fisica, Contrada Papardo, 1-98166 Messina, Italy 

^ CNR-IPCF, Viale Ferdmando Stagno d'Alcontres 37, 98158 Messina, Italy 

(Dated: July 6, 2011) 

We present a Monte Carlo simulation study of the phase behavior of two-dimensional classical 
particles repelling each other through an isotropic Gaussian potential. As in the analogous three- 
dimensional case, a reentrant-melting transition occurs upon compression for not too high tempera- 
tures, along with a spectrum of water-like anomalies in the fluid phase. However, in two dimensions 
melting is a continuous two-stage transition, with an intermediate hexatic phase which becomes in- 
creasingly more definite as pressure grows. All available evidence supports the Kosterlitz-Thouless- 
Halperin-Nelson- Young scenario for this melting transition. We expect that such a phenomenology 
can be checked in confined monolayers of charge-stabilized colloids with a softened core. 

PACS numbers: 05.20.Jj, 61.20.Ja, 64.70.D- 



In two dimensions thermal fluctuations do not allow 
the existence of a true crystalline order; in fact, only 
a quasi-long-range translational order is possible while 
bond-angular order is truly long-ranged. This opens 
the way to a two-stage melting transition through an 
intermediate "hexatic" phase with short-ranged trans- 
lational order but extended bond-angle correlations. 
In the celebrated Kosterlitz-Thouless-Halperin-Nelson- 
Young (KTHNY) theory of two-dimensional (2d) melt- 
ing [J|, the hexatic phase is promoted by the thermal 
unbinding of dislocation pairs, followed by the prolifer- 
ation of free disclinations on entering the normal fluid. 
The KTHNY theory predicts melting to be continuous. 
In two dimensions, when the energy of the dislocation 
core is sufhciently small, a first-order melting transition 
is more likely, driven by the spontaneous generation of 
grain boundaries [5|, |6|. Hexatic phases have been ob- 
served in various types of colloids ITHISI , and found also 
in some classical 



211 



J16I-|20| and quantum simulations 
Moreover, nothing prevents the hexatic phase to be just 
metastable, as observed e.g. in Ref. [2a |. 

Observing the KTHNY scenario is notoriously difficult 
because of the existence of important finite-size effects 
and long equilibration times. Also, the usually narrow 
temperature extent of the hexatic phase makes it hard to 
distinguish a two-stage melting from a single weakly first- 
order transition. Particularly severe is the situation for 
hard-core particles, where enormous samples and huge 
simulation times are required in order to_ discriminate 
between the various transition scenarios [23], while less 
demanding may be state sampling for systems of "soft" 
particles whose steric constraints are less pronounced. 

We hereby inquire into the existence of a hexatic 
phase for the 2d Gaussian-core model (GCM) pair poten- 
tial [2J], v{r) — eexp(— r^/(T^) with e > 0, which is some- 
what representative of a whole class of systems of inter- 
penetrating particles (e.g. dilute dispersions of polymer 
chains) [25|. In three dimensions, this system is known to 
exhibit reentrant melting (i.e., melting upon compression 
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FIG. 1: (Color online). Phase diagram of the 2d GCM with 
anomaly loci. Plotted in black is the melting line, with red 
dots at the computed solid-to-hexatic transition points. The 
black squares (with error bars superimposed) give the upper 
stability threshold of the solid when heated isobarically in 
steps of AT = 0.0005. The inset shows a magnified portion 
of the melting line with an adjacent (magenta) strip corre- 
sponding to the hexatic region. The dotted curves mark the 
boundary of anomaly regions (structural anomaly, red; dif- 
fusion anomaly, blue; density anomaly, black). Observe that 
the red and blue curves, which appear indistinguishable on 
the scale of the figure, depart from each other at much higher 
temperatures. 



at constant temperature) [26J as well as waterlike anoma- 



lies [27|. Except for a 30-year old canonical-ensemble 



investigation [28] with inconclusive answers, we do not 
know of any simulation study of the melting behavior of 
2d systems of particles with bounded interactions with 
a focus on the quest for a hexatic phase. Furthermore, 



it would be interesting to know about the interplay be- 
tween anomalous melting (that is, melting of the solid 
into an anomalous fluid) and the modality of decay of 
bond-angle correlations, an issue that has never been ad- 
dressed before. As discussed in more detail below, the 
melting of the 2d GCM is indeed continuous and two- 
staged, with an extremely narrow hexatic region whose 
properties comply with the predictions of the KTHNY 
theory. The complete phase diagram is plotted in Fig. 1, 
together with a number of anomaly loci in the fluid phase. 
Particles interacting through a repulsive Gaussian po- 
tential are expected to exhibit reentrant melting and a 
maximum melting temperature [29|. By examining all 
the five Bravais lattices and the honeycomb lattice, we 
first checked that the most stable state of the GCM at 
zero temperature is a triangular crystal for any pressure 
P. This gave us confidence that the triangular lattice 
provides the structure of the solid phase also for non-zero 
temperatures. We carried out isothermal-isobaric Monte 
Carlo (MC) simulations of TV-particle samples (with N 
up to 6048) in order to locate melting for a number of 
selected pressures (0.05, 0.2, 0.4, and 0.6, in reduced, 
e/cr^ units). Our method consists in running simulations 
in a sequence, starting from the cold triangular solid on 
one side of the chain and from the hot fluid on the other 
side. Then, the solid was gradually heated (the fluid 
was cooled) in temperature steps of AT = 0.0005 (in re- 
duced, e/fcs units), until we observed the abrupt melting 
(freezing) of the system. With this protocol, we found the 
same shape of the melting line as in the three-dimensional 
GCM, with a maximum melting temperature Tm of about 
0.0115 for P ^ P^ < 0.2. We plot in Fig. 1 the melt- 
ing line with three other curves which encompass re- 
gions in the fluid phase where an "anomalous" behav- 
ior occurs. On increasing the density, one first meets 
the so-called structural anomaly (which is where the ab- 
solute pair entropy reaches its maximum) [30|, followed 
by the diffusion-anomaly locus (where the self-diffusion 
coefficient [3l| attains its minimum), and by the density- 
anomaly line (where the particle-number density attains 
a local maximum). The same succession of anomaly loci 
is found in three dimensions 
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To disentangle first-order from continuous melting, we 
performed another series of runs across our earlier guess 
of the transition point, now with a 5-time larger T resolu- 
tion and also allowing for much longer equilibration times 
(10^ sweeps, that is a million MC moves per particle) 
and production runs (5 x 10^ sweeps). A typical result is 
shown in Fig. 2, where we report the average specific en- 
ergy u and particle-number density p for various system 
sizes as a function of T for P = 0.6. A continuous path 
joins the solid and fluid branches with no evidence of hys- 
teresis, which points to a smooth transition between the 
solid and fluid phases. Moreover, the energy and volume 
histograms have a simple Gaussian shape with no trace 
of bimodality within the relevant temperature range. As 
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FIG. 2: (Color online). Total energy per particle (top) 
and particle-number density (bottom) for three different sizes 
(iV = 1152, yellow; N = 2688, red; A^ = 6048, blue) for 
P — 0.6. We show results for both heating (squares) and 
cooling (dots) trajectories (one million MC sweeps of equili- 
bration plus five million sweeps of data accumulation). The 
solid and fluid branches (dotted and solid lines, respectively) 
are computed with A'^ — 1152 and much smaller statistics. 
The arrows mark the estimated transition points (see Fig. 3) . 
The small hysteresis observed upon cooling for A'' = 2688 in- 
dicates that much longer runs are needed in order that the 
solidifying system may get rid of the extra defects. 



we are going to show in the following, the intermediate 
region between the solid and the (normal) fluid can be 
qualifled as hexatic. 

We measured two different order parameters (OP), 
which are separately sensitive to the overall translational 
and orientational triangular order, with their respective 
susceptibilities and correlation functions. The transla- 
tional OP is taken to be 
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(1) 



where the sum is over the particle labels and G is any 
flrst-shell reciprocal-lattice vector of the triangular crys- 
tal. From its very definition, it follows that V't is sizeable 
only in a triangular solid that is oriented in a way con- 
sistent with the length and direction of G. Hence, ■0r is 
only measured on heating, where memory of the original 
crystal orientation is preserved as long as the system is 
large and remains solid. We anyway checked - through 
the location of the main peaks of the structure factor - 
that the orientation of the solid never changed from one 
run to the next. A sharp drop of tpx signals the melting 
of the solid into a fluid, be it hexatic or normal; concur- 



rently, the corresponding susceptibility 
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shows a distinct peak whose location is an unambiguous 
estimate of the melting transition point. At regular in- 
tervals during the simulation, we made use of the Voronoi 
construction in order to identify the nc{i) nearest neigh- 
bors (NN) of each particle i, together with the orientation 
Snn of each neighbor bond with respect to a reference 
axis. Whence, the orientational OP follows as 
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(3) 

The orientational susceptibility xe is then defined in a 
way analogous to Eq. ([2]), with ^6(ri) replacing exp(iG • 
Tj). ^g undergoes a sudden drop at the hexatic-fluid tran- 
sition, i.e., at a temperature larger than the one where 
il>T vanishes. Finally, the local bond-angular OP ^6(ri) 
enters the definition of the orientational correlation func- 
tion (OCF): 



K{r) = p 
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(4) 

where the prime over the sum excludes i = j and 
r = |R - R'|. The KTHNY theory predicts an algebraic 
f-niT) large-distance decay of the OCF in the hexatic 
phase, which should be contrasted with the exponential 
asymptotic vanishing of angular correlations in a normal 
fluid. Another prediction of the theory is ry = 1/4 at the 
hexatic-to-normal fluid transition point. 

In Fig. 3, we plot the two OPs and susceptibilities for 
P = 0.6 (an analogous behavior was observed for all the 
other pressures). We see that ipT vanishes at a slightly 
smaller temperature than iJJq, which implies that the hex- 
atic phase is confined to an extremely narrow T interval 
not wider than 0.0002-0.0003, as also witnessed by the 
maxima of the two susceptibilities occurring at slightly 
different T values. The estimated width of the hexatic 
region compares well with the temperature range of the 
bridging region between the solid and fluid branches in 
Fig. 2. While the size scaling of xe is a clear imprint of a 
second-order hexatic-to-normal fluid transition, the solid- 
to-hexatic transition might even be first-order, were this 
not in contrast to the smooth behavior of u and p. Upon 
reducing the pressure, the width of the hexatic phase 
gradually shrinks until, for P = 0.05, it becomes com- 
parable to the temperature resolution. However, even in 
this case we tend to exclude the disappearance of the hex- 
atic phase for low pressures since this would imply the 
existence of a triple point for which we currently have 
no independent evidence. Finally, it is worth mentioning 
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FIG. 3: (Color online). Order parameters and susceptibil- 
ities for P = 0.6 in the T range across the melting transi- 
tion. Upper panels: the orientational order parameter tpe and 
its susceptibility xe for three system sizes (color codes as in 
Fig. 2). Dots and squares mark data obtained by cooling and 
by heating, respectively. Roughly, the difference between the 
two estimates gives a clue about the statistical uncertainty 
associated with each data point. Lower panels: the transla- 
tional order parameter xpT and its susceptibility xt for the 
same sizes on heating. The non-zero value of ipr in the solid 
phase is actually a finite-size effect, made possible by the use 
of periodic boundary conditions in the simulations, since for 
a 2d infinite solid quasi-long-range translational order implies 
tpT ~ 0. Moreover, xt is expected to diverge in the solid 
phase of an infinite-size system. Similar considerations apply 
for the behavior of ^pa and X6 in the hexatic phase. 



the case P = 0.2 (a pressure above Pm but outside the 
density- anomaly region), where the hexatic fluid shows a 
density anomaly while the normal fluid does not - in no 
other way could the density branch of the normal fluid 
have hooked on a solid branch that lies at a lower density 
level. 

A more direct evidence of the hexatic phase emerges 
from the large-distance behavior of the OCF. We plot 
this function in Fig. 4 at various temperatures across the 
hexatic phase for P — 0.6. It appears that the OCF de- 
cays algebraically in a T region of limited extent, which 
roughly corresponds to the middle of the bridging region 
in Fig. 2. Moreover, the decay exponent in this hexatic 
region is smaller than 1/4, becoming larger only on pass- 
ing to the normal fluid. 

We finally checked a further KTHNY prediction con- 
cerning the behavior of a 2d triangular solid which is 
about to melt into a hexatic phase. The elastic insta- 
bility that signals the onset of dissociation of dislocation 
pairs, preluding to the stabilization of the hexatic phase. 
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FIG. 4: (Color online). Orietational correlation function 
hs(r) at selected temperatures across the hexatic region for 
P = 0.6. We plot hs{r) on heating for two sizes, N = 2688 
(red) and N = 6048 (blue). Top: log-log plot; bottom: log- 
lin plot. Upon increasing T from 0.0054 to 0.0059 there is 
a qualitative change in the large-distance behavior of he{r), 
from constant (solid) to power-law decay (hexatic fluid), up 
to exponential decay (normal fluid). Note that, consistently 
with the KTHNY theory, the decay exponent 77 is less than 
1/4 (i.e., the slope of the dotted curve) in the hexatic phase. 



is heralded by the value of 



K = 



4a^ fi{l-j. + A) 



(5) 



becoming equal to 167r [J, |3J]. In Eq. ([S]), A and /x are 
the Lame coefficients (as renormalized by the thermal 



fluctuations) while a = y 2/{\/3p) is the lattice param- 
eter. A and /i are respectively given by C12 + P and 
C44 — P, in terms of the elastic constants C12 and C44 
which can be computed as canonical-ensemble averages 
from virial-like formulae 3J|. We found an impressive 
confirmation of the theory for P = 0.6 while the thresh- 
old value of K/{1Qtt) (before its drop to zero) turned out 
to be a bit larger than one (1.1-1.2) for the other inves- 
tigated pressures. This further indicates that the over- 
all KTHNY picture deteriorates with reducing pressure, 
probably because the formation energy of a dislocation 
becomes smaller and smaller with increasing average in- 
terparticle distances. 

In conclusion, we have provided the first unambiguous 
evidence of the occurrence of two-stage continuous reen- 
trant melting via a hexatic phase in the 2d Gaussian- 
core model, taken as prototypical of the phase behavior 
of bounded model potentials. We have validated a num- 
ber of KTHNY predictions, though larger samples and 



more statistics will be necessary in order to ascertain 
the real nature of melting at low densities. The present 
discovery of reentrant-hexatic behavior in the GCM is 
relevant for many soft-matter systems. For instance, 
one can engineer colloidal particles interacting through 
a temperature-modulated softened repulsion, which will 
likely exhibit GCM-like reentrant melting in a range of 
packing fractions well below the density at which hard- 
core crystallization occurs (see [35j for a 3d realization 
of this scenario). Such systems would be natural can- 
didates where to detect (e.g. by video microscopy [36]) 
a reentrant-hexatic phenomenon of the kind illustrated 
here. 
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